Fast Inverse Square Root

2016-02-02

前言

Fast Inverse Square Root 即平方根倒数速算法，是一种快速计算$x^{\frac{1}{2}}$ 的方法。在公布的《雷神之锤3》里被发现，写的很强大，原作者不详，但是与卡神相关。遇到这种神级代码，就想记录下来。_(:зゝ∠)_

Code

先来看看这段代码，其中的注释也是原汁原味的。

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;
    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking（对浮点数的邪恶位级hack）
    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?（这他妈的是怎么回事？）
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration （第一次牛顿迭代）
//      y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed（第二次迭代，可以删除）
    return y;
}

魔术数字

平方根倒数速算法在速度上的优势源自将浮点数转化为长整型以作整数看待，并用特定常数0x5f3759df与之相减。然而代码中出现的这个0x5f3759df十分诡异，着实让人摸不着头脑。因此它跟出现在其他代码中让人难以理解的常数一样，被称为魔术数字。

将浮点数当作整数先位移后减法，所得的浮点数结果即是对输入数字的平方根倒数的粗略估计值，而后再进行一次牛顿迭代法，以使之更精确后，代码即执行完毕。

为何是0x5f3759df?

Chris Lomont曾做了个研究：他推算出了一个函数以讨论此速算法的误差，并找出了使误差最小的最佳R值0x5f37642f（与代码中使用的0x5f3759df相当接近），但以之代入算法计算并进行一次牛顿迭代后，所得近似值之精度仍略低于代入0x5f3759df的结果；因此Lomont将目标改为查找在进行1－2次牛顿迭代后能得到最大精度的R值，在暴力搜索后得出最优R值为0x5f375a86，以此值代入算法并进行牛顿迭代，所得的结果都比代入原始值（0x5f3759df）更精确。取自weki平方根倒数速算法。

牛顿迭代法

y  = y * ( threehalfs - ( x2 * y * y ) );

这行代码来源于牛顿迭代法（或叫牛顿方法）。具体推到如下：

设：
$$
y = \frac{1}{\sqrt{x}}
$$

将其化作以$y$做自变量。

$$
f(y) = \frac{1}{y^2} - x
$$

则有：

$$
f^{‘}(y) = \frac{-2}{y^{3}}
$$

所以由 牛顿迭代公式 得：

$$
y_{n+1} = y_{n} - \frac{f(y_{n})}{f^{‘}(y_{n})}\\
= \frac{y_{n}(3-xy_{n}^{3})}{2}\\
= y_{n}(1.5 - \frac{xy_{n}^{3}}{2})
$$

即对应代码中的该行：

y  = y * ( threehalfs - ( x2 * y * y ) );